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ABSTRACT 


^  An  approach  to  modelling  and  residual  analysis  of  nonlinear 
autoregressive  time  series  in  exponential  variables  is  presented;  the 
approach  is  illustrated  by  analysis  of  a  long  series  of  wind  velocity 
data  which  has  first  been  detrended  and  then  transformed  into  a  stationary 
series  with  an  exponential  marginal  distribution.  The  stationary  series  is 
modelled  with  a  newly  developed  type  of  second  order  autoregressive  process 
with  random  coefficients,  called  the  NEAR( 2 )  model;  it  has  a  second  order 
autoregressive  correlation  structure  but  is  nonlinear  because  its 
coefficients  are  random.  The  exponential  distributional  assumptions 
involved  in  this  model  highlight  a  very  broad  four  parameter  structure 
which  combines  five  exponential  random  variables  into  a  sixth  exponential 
random  variable;  other  applications  of  this  structure  are  briefly 
consioared.  Dependency  in  the  NEAR( 2  >  process  not  accounted  for  by 
standard  autocorrelations  is  explored  by  developing  a  residual  analysis  for 
ti:<te  series  having  autoregressive  correlation  structure;  this  involves 
defining  linear  uncorrelated  residuals  which  are  dependent,  and  then 
assessing  this  higher  order  dependence  by  standard  time  series 
computations.  The>  application  of  this  residual  analysis  to  the  wind 
velocity  data  illustrates  both  the  utility  and  difficulty  of  nonlinear  time 
series  modelling.  “V?. 


1 


i.  imtboouction 

There  are  several  aspects  of  many  observed  univariate  time  series  which 
are  not  satisfactorily  accounted  for  in  standard  time  series  analysis*  they 
include  nonGaussian  marginal  distributions,  dependence  not  accounted  for  by 
second  order  moments  (autocorrelations)  and  directionality  in  the  time 
series.  Quite  often  a  Gaussian  distribution  will  be  inappropriate  because 
the  variable  being  modelled  has  a  positive  and  highly  skewed  distribution, 
e.g.  wind  speeds,  the  service  times  in  a  queue,  or  the  daily  flows  of  a 
river.  Many  particular  such  distributions  can  be  envisaged  and  time  series 
models  have  been  constructed  for  them.  Examples  are  Gamma  distributions 
(Gaver  and  Lewis,  1980;  Lewis,  1981;  McKenzie,  1982,  Lawrance,  1982)  and 
mixed  exponential  distributions  (Gaver  and  Lewis,  1980;  Lawrance,  1980a; 
Lawrance  and  Lewis,  1982). 

However  the  simplest,  most  widely  used  and  most  analytically  tractable 
of  these  distribution  models  is  the  exponential  distribution.  Like 
Gaussian  random  variables,  exponentially  distributed  random  variables  enjoy 
many  special  properties;  also  they  can  be  mildly  transformed  quite  easily 
into  distributions  which  are  either  more  skewed  or  less  skewed  than  the 
exponential.  The  Weibull  distribution  is  an  example,  being  just  a  power 
transformation  of  an  exponentially  distributed  random  variable.  Thus  the 
approach  here,  following  earlier  work  (Gaver  and  Lewis,  1980;  Lawrance  and 
Lewis,  1080 ,  1981)  is  to  regard  the  exponential  variables  as  canonical  and 
tc  develop  their  use  in  time  series  modelling. 

It  should  also  be  noted  that  time  series  of  (marginally)  uniformly 
distributed  random  variables  can  be  obtained  by  exponential  transformations 
of  time  series  in  exponentially  distributed  variables;  such  uniform 
processes  could  then  be  used  to  generate  time  series  with  other  desired 


marginal  distributions. 


The  work  cited  previously  has  concentrated  for  the  most  part  on  first 
order,  nonGaussian  autoregressive  models,  both  of  the  standard  type 
(constant  coefficient,  additive,  linear  combinations)  and  a  random 
coefficient  type  Introduced  by  the  authors.  The  extension  of  the  models  to 
higher  order  autoregression  is  clearly  necessary  to  attain  flexibility  in 
modelling  correlation  and  dependency  structure  of  the  processes,  but  these 
extensions  are  in  no  way  as  immediate  as  in  the  standard  linear  Gaussian 
case.  A  simple  mixing  device  cam  be  used  (Jacobs  and  Lewis,  1983)  but  the 
range  of  correlations  attained  is  much  narrower  tham  the  range  attained  in 
the  standard  linear,  second  order  autoregressive  structure.  A  broader 
extension,  called  the  EAR( 2 )  model,  yam  obtained  in  the  exponential  case  by 
Lawrance  and  Lewis  (1980),  but  its  innovation  variable  has  a  zero  component 
which  gives  runs  in  the  process;  this  will  often  be  hard  to  justify. 

A  major  part  of  the  present  work  consists  of  obtaining  a  very  broad  and 
rich  extension  of  the  NEAR(l)  model  (Lawrance  and  Lewis,  1981)  to  a  second 
order  autoregressive  process;  it  includes  the  EAR(  2 )  model  but  does  not 
generally  have  a  zero  component.  This  NEAR(2)  model  was  proposed  in 
Lawrance  (1980b),  later  reviewed  in  Raftery  (1981),  but  the  necessary 
analysis  of  its  innovation  structure  was  not  given.  Here  the  innovation 
random  variable  for  the  NEAR(  2 )  process  is  proved  to  exist  without 
unnatural  boundaries  on  its  (four)  parameter  region;  explicit  construction 
is  given  for  the  innovation  random  variable. 

The  richness  of  the  four  parameter  NEAR( 2 )  model,  and  the  fact  that  an 
infinite  number  of  cases  of  the  model  with  identical  correlation  structure 
are  available,  forces  consideration  of  higher  order  aspects  of  dependence. 
The  analysis  of  the  higher  order  aspects  of  exponential  time  series  is  at  a 
fairly  early  stage  and  is  as  follows.  First  it  will  be  shown  that  the 
autocorrelations  p(t),  I  -  0,*1,±2,...  for  the  NEAR(2)  process  satisfy  the 
Yule-Walker  equations  with  constants  ai  and  a?  Which  are  functions  of  the 


four  parameters  of  the  model.  This  fol lows  immediately  from  the  fact  that 
Xfx  is  a  random  ceofficient,  linear  additive  combination  of  Xn-^, Xn_2  and 
the  innovation  random  variable  en.  Secondly,  it  can  be  shown  (Lewis  and 
Iawrance,  1984)  that  the  residuals  Xn  -  ajXn-i  -  a2Xn_2»  Which  are  the  usual 
residuals  for  second  order  constant  coefici.ent,  linear  additive 

autoregressive  processes  are  uncorrelated. 

Thus,  although  the  standard  analysis  of  time  series  stops  with 
uncorrelated  residuals,  i.e.  a  flat  spectrum  for  the  residuals,  such 
residuals  can  also  be  used  to  good  effect  to  investigate  higher  order 
aspects  of  dependence  in  the  NEAR( 2 )  model,  in  fact,  if  the  autoregression 
is  not  of  the  standard  type  ( constant  coefficient,  additive,  linear 
combinations)  the  (uncorrelated)  residual  will  be  dependent.  One  aspect  of 
this  is  that  the  squared  residuals  will  have  non-zero  autocorrelations,  and 
another  is  that  the  crosscorrelations  of  residuals  and  squared  residuals 
will  be  non- zero;  both  sets  of  correlations  are  theoretically  zero  when  a 
standard  second  order  autoregressive  model  is  appropriate. 

This  residual  analysis  will  be  illustrated  by  some  theoretical 
calculations  for  the  NEAR(  2 )  model  and  by  a  brief  application  to  a  long 
series  of  detrended  and  transformed  wind  velocity  data. 


Our  aim  in  this  section  is  to  give  in  outline  the  ideas  leading  to  the 
tine  series  model  of  main  concern  in  this  paper,  and  called  NEAR(2), 
following  the  earlier  terminology  NEAR(l)  in  Lawrance  and  Lewis  (1981). 
The  NEAR( 2 )  model  has  four  parameters,  and  incorporates  and  broadens  the 
earlier  two  parameter  EAR( 2 )  model  (Lawrance  and  Lewis,  1980).  The  NEAR( 2 ) 
model  will  be  exponential  in  marginal  distribution,  have  second  order 
autoregressive  Markov  dependence,  and  have  autocorrelations  satisfying 
second  order  difference  equations  of  the  familiar  Yule-Walker  type.  In 
addition  it  will  have  dependence  beyond  autocorrelation,  and  will  not  be 
reversible  in  time.  It  is  not  linear  in  the  standard  sense,  having  random 
coefficient,  linear  additive  autoregressive  structure,  but  neither  is  it 
nonlinear  in  the  standard  sense  of  incorporating  powers  or  products  of 
lagged  variables.  Also,  it  differs  from  the  random  coefficient  models 
considered  by  Nicholls  and  Quinn  ( 1982 )  in  that  the  marginal  distribution 
is  specified.  The  view  taken  here  is  that  the  marginal  distribution  is  the 
easiest  aspect  of  data  to  look  at  and  should  be  the  starting  point  for 
modelling. 

Writing  (Xn)  for  the  time  series  variables,  and  (En)  for  an  i.i.d. 

exponential  innovation  sequence  of  unit  mean,  the  two  parameter  NEAR(l) 

model,  as  previously  defined,  is  given  by 

f*n-l  w.p.  a  En  w.p.  p 

X„  -  +  ,  (2.1) 

.  0  w.p.  1-a  tbEn  w.p.  1-p 

wit’i  b-(l-a)/3  and  p-(  l-0)/{l-(  l-a)0) .  The  parameter  region  is,  in  general, 
O<a,0<l,  l.  The  case  0-1,  o«x<l  is  rather  special,  and  has  been  called 

the  TRAR(1)  model,  and  when  a-1,  0</3<l,  the  earlier  EAR(l)  model  is 
recovered.  Except  for  this  KAR( 1 )  case,  the  NEAR(l)  model  does  not  allow 
zero  innovations  (Caver  and  Lewis,  1980)  and  so  is  more  statistically 
acceptable.  The  zero  innovation  in  the  EAR(  1 )  case  implies  that  Kn*(fX.n-i 


and  thus  0  can  be  determined  exactly  from  runs  down  jn  the  sample  path  of 


the  process. 

in  general  the  i.i.d.  innovations  in  the  NEAR(  1 )  process  are  formed  as 
the  probabilistic  mixture  of  two  exponentials,  and  are  thus  easily 
simulated. 

The  NEAR(2)  model  is  a  direct  generalization  of  (2.1)  and  takes  the 
form 

‘01xn-l  w-P-  ®1 

Xn  =  02xn-2  W-P-  a2  +  en  (2.2) 

0  w.p.  l-ai-o<2 

with  parameter  region  a^O ,a,2*0, a1+a2<l , O<0i , 02<1  J  («nJ  is  an  appropriately 
chosen  innovation  sequence.  Many  special  cases  can  arise  when  the  above 
restrictions  include  some  of  the  equalities  and,  for  the  purposes  of  a 
general  development,  it  is  best  to  regard  the  inequalities  as  strict.  Given 
that  (Xn)  is  required  to  have  an  exponential  marginal  distribution,  the 
main  question  concerns  whether  there  is  a  valid  probability  distribution 
for  en.  The  Theorem  proved  in  Section  2.3  will  show  that  this  is  the  case, 
and  that  the  distribution,  when  the  inequalities  on  a^,a2  and  0^,02  the 
parameter  region  are  strict,  takes  the  form 

En  w.p.  I-P2-P3 

en  -  b2En  w.p.  P2  •  (2.3) 

b3En  w.p.  p3 

a  probabilistic  mixture  of  three  exponentials  with  parameters  given  in 
Section  2.3.  To  establish  this  result  a  fairly  detailed  analysis  of  a 
derived  moment  generating  function  is  required.  This  is  necessary  since  a 
direct  moir.ant  generating  function  solution  of  (2.2)  for  en  does  not 
establish  chat  en  has  a  proper  distribution;  all  that  is  shown  is  that  the 
solution  is  a  possibly- improper  mixture  of  three  exponentials. 


3.  VALIDITY  or  THE  HEAR<2  )  MODEL 

In  this  section  we  prove  the  following 
THEOREM.  Let  (En)  be  am  i.i.d.  sequence  of  unit  mean  exponential  random 
variables.  Then  if  the  four  parameters  oi,O2,0i»02  satisfy  ai>0,O2>0, 
al+a2<1»o<01'02<1'  the  relationship 
OlXn-l  W*P*  «1 

Xf]  *  02^1-2  w . p .  Gig  ^  «n'  n-O ril^ t2f .  ( . f  (3.1) 

.0  w.p.  l-oti-0(2 

where 

En  w.p.  I-P2-P3 

«n  "  b2En  w.p.  P2  ,  ( 3 • 2 ) 

b3En  w.p.  p3 


defines  a  stationary  sequence  of  (marginally)  exponentially  distributed 
random  variables  with  mean  one.  Here 


P2  “  ( ( a101+a202 )b2  ~  ( «i+a2 )01^2 )/(( b2~b3 )( l-b3 ) )  , 
p3  *  { ( <*i+<*2  )0102  ~  (ai<3i+«2<32)b3}/((b2~b3)(1-b3))  * 


(3.3) 


(3.4) 


o  <  b3  -  (s-(  s*-4r  )l/z)/2  <  b2  -  {s+ts^r)1/*}/^  <  1 


where 


s  *•  ( l-ai  )0i  +  ( l-a2  )^2 
r  »  ( l-ax~a2 )Pl^2 • 


(3.5) 


(3.6) 


(3.7) 


proof.  For  the  NEAR(2)  model  specified  by  (3.1)-(3.7),  let  an<s  $«(t) 

be  the  moment  generating  functions  of  the  (Xn)  amd  (cn)  sequences;  then  if 


stationarity  of  the  (Xn)  series  is  assumed, 

$X<t)  -  ^e(t){ai^x(0xt)  +  a2^x(02t)  +  ( 1-01-02)) 


(3.8) 


Assuming  an  exponential  marginal  distribution  of  unit  mean  for  (Xn),  then 
the  independent  distribution  of  (en)  hats  moment  generating  function, 
po£3ibly  not  proper,  given  by 

A  _  _  _ _ Ll±Pl*lU:±PZ*l _  ,3  g) 

'  (l+t)t(l-o1-o2)0i02t*  +  {(l-o1)01  +  (l-a2)02)t  +  1]  '  K  ’ 

It  is  convenient  to  establish  right  away  that  the  quadratic  term  in  the 
denominator  of  (3.9)  has  real  distinct  and  positive  roots,  bi  and  b2;  this 
eliminates  any  subsequent  need  to  invert  such  a  term  ats  a  whole.  The 


required  condition  for  real  distinct  roots  is  that 
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[(1-Oi)0i  +  (l-o<2>02]*  -  4(l-ai-a2)0i02 
be  positive.*  this  is  so  from  its  equality  to  the  expression 

t(l-Oi)0i  -  (1-02 )02l 2  +  4aia20i02 

which  is  clearly  positive;  the  positivity  o£  the  roots  bi  and  b2  is  obvious 
from  (3.9)  since  their  product  and  sum  given  in  (3.11)  and  (3.12)  below  are 
both  positve. 

With  bi  and  b2  positive,  a  partial  fraction  expansion  of  (3.2)  can  be 
written  in  the  suggestive  form 

*«(t>  "  <  1-P2-P3 )  "t  +  P2  Ii£it  +  P3  1+SJt  •  (3-10) 

Comparisons  between  (3.9)  and  (3.10)  then  show  that  b3,b3  and  P2»Pa  may  be 

obtained  in  terms  of  0i,02  and  01,02  by  solving  the  equations 


b2  +  b3  =  (1-Oi)0i  +  (1-02)02*  (3.11) 

b3b3  =  ( 1-01-02 )0i02  *  (3.12) 

(l-b2)P2  +  (  l~b3  )P3  »  <*101  +  <*202*  (3.13) 

b3( l-b2 )P2  +  b2( l-b3 )P3  =  ( 01+02  )0i02  •  (3.14) 


A  difficulty  with  this  apparently  straightforward  solution  is  that  the 
inversion  of  (3.9)  or  (3.10)  could  lead  to  a  function  which  is  not  a 
probability  density,  or  it  could  yield  a  probability  density  but  not  one 
which  is  a  probabilistic  mixture  of  three  exponentials.  In  fact,  neither  of 
these  possibilities  is  the  case,  as  will  be  shown  by  establishing  that  p2 
and  p3  are  positive  and  subject  to  the  condition  P2+P3<1*  and  hence  can 
represent  probabilities. 

Explicit  expressions  for  P2  and  P3  can  be  obtained  from  (3.13)  and 
(3.11)  and  are  given  at  (3.3)  and  (3.4).  Prom  now  on  it  will  be  assumed,  in 
accordance  with  the  theorem,  that  b2  is  the  larger  of  b2  and  b3,  these  being 
obtained  by  solving  the  quadratic  pair  (3.11)  and  (3.12).  To  establish  that 
P2+P3<1,  we  have,  by  adding  (3.3)  and  (3.4), 

p-+p,  -  >  I_<_«lt«2 )£l?2  .  (3.15) 

P2  P3  (l-b2)<l-b3)  11 

Multiplying  out  (l-b2)(l-b3)  in  the  denominator  and  using  (3.11)  and  (3.12) 


gives,  after  some  rearrangement. 


P2+P3  -  1  - 


(3.16) 


(  1-/3!  )( l-/32  )  +  Oi/3!(  1-/32  >  +  a2/32(  1-Pi  > 

The  algebraic  expression  here  is  clearly  positive  and  less  than  one,  from 
which  it  follows  that  Pi+P2<l. 

The  positivity  of  p2  and  P3  will  now  be  proved  by  showing  that  the 
numerators  and  denominators  of  (3.3)  and  (3.4)  are  positive.  For  the 
denominators,  this  requires  that  0<b2,b3<l  which  will  be  verified  by 
showing  that  0<b2b3<l  and  0<( l-b2 )( l-b3 )<1.  The  first  of  these  latter  two 
inequalities  is  obvious  from  (3.12);  for  the  second  consider  the 
expressions 


(l-b2)(l-b3)  =  1  -  (b2+b3)  +  b2b3 

-  ( ai0i+l-0i )(  a2/32+l-/32  )  -  ( ax0x )( a2&2  )  (3*17) 

after  using  (3.11)  and  (3.12),  and  then 

1  -  (l-b2)(l-b3)  *  b2  +  b3  -  b2b3 

-  (  l-ax  )0X(  l-/32  )  +  (  l-a2  )/32(  1-/3X  )  +  0X02  •  (3.18) 

The  right  hand  sides  of  both  (3.17)  and  (3.18)  are  obviously  positive.  This 
concludes  the  proof  that  0<b2,b3<l  and  hence  that  the  denominators  of  p2  and 
P3  are  positive. 

For  the  numerators  of  p2  and  P3  to  be  positive  (3.3)  and  (3.4)  indicate 

that  b  =  (ax+a2)px02/(axpx+a202)  must  satisfy  the  inequalities 

b3  <  b  <  b2.  (3.19) 

At  this  last  stage,  explicit  expressions  for  b2  and  b3  must  be  used,  and 

from  (3.11)  and  (3.12)  are  given,  after  writing 

s  =  (l-ax)/3i  +  ( l-a2  )/32  and  r  «  ( l-ai~a2  )px02  > 
by  (3.20) 

b2  *  (s  +  (82-4r)x/2}/2  and  b3  *  {s  -  ( s2-4r  )A/2}/2 

Then  (3.19)  is  equivalent  to 

-  (s2-4r)1/2  <  s  -  2b  <  (s2-4r)1/2 


or 

s2-4r  >  (s-2b)* 

or 

sb  -  b2  -  r  >  0  * 

(3.21) 

After  some  algebraic  rearrangement  the  left  hand  side  of  (3.21)  becomes 

ala20102<  01“02  )*/(  °ldl+<*2^2  )*  (3.22) 
Which  is  again  clearly  strictly  positive,  as  was  to  be  proved. 

This  concludes  the  proof  that  P2  and  P3  are  both  positive  and  subject  to 
p2+p3<l;  hence  l-p2-P3 ,  Pz>  and  P3  can  all  be  regarded  as  probabilities. 
Thus  «n  has  a  proper  probability  distribution  which  can  be  generated  as  the 
( l~P2~P3'P2'P3 )  mixture  of  three  exponentials  of  means  1,  b2  and  b3 
respectively;  further,  both  b2  and  b3  are  less  than  unity  and  b2  *  b3 . 

In  special  cases  there  are  valid  and  simpler  results  for  the 
distribution  of  en.  For  instance,  when  01*02*1 »  <n  bas  a  simple  exponential 
distribution  of  mean  (1-01-02).  When  /3i»02^1  the  innovation  has  a  mixed 
exponential  distribution  of  the  NEAR( 1 )  form  given  in  (2.1)  with  o*oi+o2 • 
When  02*1'  P2+P3*l  and  «n  is  the  mixture  of  two  exponentials  with  means  b2 
and  b3;  this  case  is  used  in  some  of  the  calculations  of  Section  9. 
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4.  OTHER  OSES  OP  THE  MEAR(2) 


The  NEAR(  2 )  process  was  established  by  showing  that  (3.2)  was  a  valid 
innovation  distribution  for  the  relation  (3.1)  to  give  a  process  with 
marginal  exponential  distributions.  The  distributional  assumptions  implied 
by  this  result  can  also  be  taken  out  of  the  time  series  context  in  which 
they  were  derived  and  viewed  generally  as  a  way  to  combine  a  pair  of 
(possibly  dependent)  unit  exponential  variables  ( Lx . L2 )  with  an  independent 
triple  of  possibly  dependent,  unit  exponential  variables  (Mx.M2.M3)  so  as 
to  yield  a  further  unit  exponential  variable.  Specifically,  with 
(al'a2'01'02)  and  (b2.b3.p2.P3)  as  previously  related  by  (3.3)  -  (3.7),  the 
Theorem  has  established  that 


Sll.1  w.p.  ax  Mx  w.p.  I-P2-P3 

02^*2  W*P*  a2  +  ^2**2  W.p.  P2 

0  w.p.  l-ax~a2  U>3**3  w.p.  P3 


(4.1) 


has  a  unit  exponential  distribution. 

First  of  all,  the  idea  of  "switching"  will  be  illustrated;  in  the 
NEAR( 2 )  context,  this  suggests  talcing  (Mx.M2.M3)  as  (Xn-x.Xn_2.xn-3 )  and 
( Lx . I«2 )  as  (En,En).  Then  (4.1)  gives  the  time  series  model 


xn-l  w.p.  I-P2-P3  [0lEn  w.p.  ax 

xn  *  b2xn-2  w.p.  p2  +  02En  w.p.  a2 

b3Xn_3  w.p.  p3  0  w.p.  l-ai-a2 


(4.2) 


This  is  a  third  order  autoregression,  actually  a  case  of  the  EAR( 3 )  model 
cited  in  Lawrance  and  Lewis  (1980);  note,  however,  that  this  third  order 
autoregressive  exponential  process  allows  zero  innovations.  another, 
better  behaved  higher  order  exponential  model  -  in  fact  a  p-th  order  model  - 
is  obtained  by  the  following  application  of  the  result  (4.1)  in  its  original 
form  (3.1).  Let  the  indices  l,2,...,p  be  partitioned  into  two  non-empty 
sets  lx  and  I 2  of  size  tx  and  t2  respectively.  Then  in  the  model 


Xn  - 


f01,xn-l  W-P-  «l’ 


0p'Xn-p  w.p.  ap- 
l  0  w.p.  1-ax*... -Op' 


n»0,*l,*2, .... 


(4.3) 


!  4.VO  .  *  a*  a*  t  *  .  *  a*  ./a*  a*  a*  1.  •  a*  .  *  a*  a'  a*  .  *  a'  *  \  *■’-  •*  .V  «*.  ** 
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let  ieli ;  0±'^02’  i«*2>  jEj  «i*“«i  iii2ai  *“2-  Th«n  if 

al+a2<1<  o<Oi,02<1»  the  distribution  of  En  is  given  by  the  Theorem.  Thus  we 

have  a  pth  order  exponential  autoregressive  process  with  four  parameters . 

However,  while  this  may  seem  satisfying  it  is  not  clear  that  four  parameters 

would  be  sufficient  to  characterize  the  sample  path  behaviour  of  an 

exponential  process  with  very  high  order  dependence. 

Another  use  of  (4.1)  is  to  allow  Lx  and  L2  to  both  be  Xn-x,  and  so 
obtain  a  four  parameter  first  order  model  of  the  form 

PlXn-l  W*P*  “1  En  W-P*  1~P2“P3 

Xn  -  02xn-l  W*P*  a2  +  b2En  W‘P*  P2  •  (4.4) 

0  w.p.  l-ot2 "a2  b3En  w.p.  P3 

Four  parameters  may  seem  excessive  for  a  first  order  autoregressive  process 
but  simulations  show  a  wide  range  of  behaviour  in  sample  paths  with 
different  choices  of  parameters.  Equation  (4.4)  in  turn  suggests  a  first 
order  model  allowing  negative  dependence.  This  is  obtained  by  replacing 
the  variable  Xn_x  in  (4.4)  which  is  multiplied  by  02  by  the  antithetic 
transformation  of  Xn_x,  which  is  log(l-exp( -Xn_i )) .  Two  parameter  versions 
of  these  two  first  order  models  could  be  obtained,  for  example,  by  taking 

alma2 *  0im&2 • 

A  third  type  of  use  of  the  construction  is  to  give  mixed  autoregressive 
moving  average  models;  for  this,  (Lx,  1*2)  is  (Xn_x»Xn-2)  48  in  (3.1),  but 
( Mx » M2 » M3 )  are  chosen  to  be  ( En,En+x*En+l )  for  a  second  order  moving 
average  component,  or  as  (En,En+X’En+2 )  f°r  a  third  order  moving  average 
component;  these  forward  running  indices  of  the  innovation  sequence  are 
necessary  for  the  required  independence  in  the  construction. 

another  use  of  this  construction  is  to  obtain  an  explicit  mixed  first 
order  autoregressive  moving  average  exponential  process,  which  could  be 
contrasted  with  the  implicit  model  given  in  Jacobs  and  Lewis  (1977).  Thus 
in  the  basic  structure  (4.1)  replace  Lx  by  Xn-x*  I»2  En-1  and  Mx.M3.M3 
each  by  En;  although  Xn_x  and  En_x  are  dependent,  they  are  independent  of  En 
as  required  by  the  construction. 


Out  of  the  tune  series  context,  the  construction  suggests  ways  to 
obtain  Multivariate  exponential  distributions,  rather  as  in  Lawrance  and 
Lewis  (1983). 

further  possibilities  are  numerous,  but  it  is  not  the  intention  here  to 
exhaustive ly  list  them,  or  to  derive  the  details  of  those  cited.  Analysis 
in  the  following  sections  will  deal  with  the  basic  NEAR(2)  model. 


S.  AUTOCORRELATION  STRUCTURE  OF  THE  NEAR(2)  PROCESS 

In  this  section  we  show  that  the  autocorrelations  p( I )»corr(Xn,Xn-f ), 
8*0,  ±l,  ±2, . . .  of  the  NEAR(  2 )  process  satisfy  AR(2)  Yule-Walker  type 
difference  equations;  thus  the  second  order  dependency  of  the  process  is 
indistinguishable  from  that  of  the  standard  autoregressive  model,  AR(2). 
To  show  this,  it  is  convenient  to  write  the  equation  (3.1)  as  a  random 
coefficient  additive  combination  of  Xn_x,  xn_2  and  En.  Thus  we  have  the 
NEAR(2)  process  in  its  random  coefficient,  linear,  additive  form  as 

Xn  *  01Kn ' xn-l  +  /32Kn"Xn_2  +  LnEn  n-0,  ±1,  *2, . . . ,  (5.1) 

where 

1  w.p.  l-p2-P3 

Ln  -  b2  w.p.  p2  ,  n-0,±l,*2, . . .,  (5.2) 

>3  w.p.  p3 

[(1/0)  w.p.  ax 

(Kn’,Kn")  =  (0,1)  w.p.  a2  ;  n-0,tl,*2, . . .  (5.3) 

(0,0)  w.p.  l-ax~a2 

the  i.i.d.  sequences  (Ln)  and  (Kn‘  ,Kn" )  are  assumed  to  be  mutually 
independent  and  Independent  of  the  independent  exponential  sequence  (En); 
the  En  *  s  are  assumed  to  have  unit  mean,  as  then  do  the  Xn’s  by 
construction . 

Now  E(Kn'  )-ax  and  E(Kn”)=a2,  80  that  E(  )-l-ai0x-a202 •  Then 

multiplying  Xn  in  (5.1)  by  Xn_f  we  have,  for  8  >  1, 

E(XnXn-8>  *  «10lE(Xn-iXn-|)  +  a202E(Xn-2Xn_4)  +  E(  L„  )E( En )E( X„_ , ) 

=  ai0iE<  Xn_  jXn_  j  )  +  a202E(  xn-2Xn- 8  )  +  1_a101~a202 » 

so  that 

E(XnXn_j)-l  -  al01{E(Xn_1Xn_J)-l}  +  a202{E(Xn_2Xn-|)-l) 
and  thus 

p(l)  «  ax  +  a2p(i),  p(2)  =  axp(l)  +  a2,  (5.4) 

0  P(  8  )  -  axp(  8-1 )  +  a2p(  8-2  ),  8-3,4,..., 

where  ax»ax0x  and  a2-ct2/32.  The  equations  (5.4)  are  the  same  as  Yule-Walker 

equations  for  the  standard  AR( 2 )  process .  The  conditions  for  a  solution  to 

exist,  namely  ax+a2<l,  ax-a2>-l,  a2>-l  are  clearly  satisfied  when  the 

conditions  on  <*x , a2  *  Pi  >  02  given  in  the  Theorem  of  Section  3  hold. 


Note  from  (5.4)  the  explicit  results 


p(l)  =  ai/(l-a2)  and  p(2)  -  axp(l)  +  a2.  (5.5) 

and  hence,  since  0<ax,a2<i,  the  restriction  of  the  autocorrelations  p(t)  to 
positive  values.  The  possible  region  of  (p(l),p(2))  values  is  bounded 
below  by  p(2)=(p(l))2  and  otherwise  bounded  by  p(l)>0  and  p(2)<l. 
Broadening  of  the  model  to  negative  dependency  may  be  achieved  using 
antithetic  ideas,  or  the  bivariate  scheme  given  in  Gaver  and  Lewis  (1980), 
but  is  not  pursued  here. 

Note  too  that  the  parameters  in  (5.4)  enter  only  as  products  ax^-aidi  and 
a2^a2/32.  Thus  for  small  enough  ax  and  a2,  values  of  Pi  and  02  greater  than 
unity  could  be  allowed,  and  (5.4)  would  still  have  a  stable  solution. 
However,  the  sequence  <sn  in  the  defining  equation  (3.2)  may  not  exist;  it 
has  not  been  determined  whether  (0i<1,0241)  is  a  necessary  condition  for 
this  existence. 

specifying  allowable  values  of  p(l)  and  p(2),  as  may  be  done  in  an 

initial  second  order  analysis  of  data,  leaves  two  parameters  to  be 

specified  in  the  model,  say  ax  and  a2,  which  could  produce  very  different 

sample  path  behaviour  in  the  time  series.  It  is  important  to  notice  that 

this  specification  of  p(l)  and  p(2)  further  constrains  the  range  of 

possible  ax  and  a2  values.  Recalling  that  p(l)  and  p(2)  fix  ax=«i/3i  and  *2 

=  a202,  as  well  as  that  ai+a2<l,  it  is  easily  shown  that  we  must  have 

al  <  al  and  a2  <  a2  (5.6) 

which  implies  that  ax+a2<ax+a2<l .  Thus  ax  and  a2  are  forced  to  lie  in  a 

triangular  subregion  of  the  triangular  (ax,a2)  region  which  is  bounded 

bel^w  by  a2,  bounded  on  the  left  by  ax,  and  bounded  above  by  the  line 

al+a25*^ •  These  results  are  useful,  and  will  be  employed  in  an  exploratory 


analysis  of  the  wind  velocity  data  in  Section  9. 


6.  Ml  ANALYSIS  OP  A  LONG  SERIES  GT  NXMD  VELOCITY  DATA 
6.1  Discussion  of  the  Data 

Lewis  and  Hugus  (1982)  have  given  an  analysis  of  a  set  of  43,800 
3-hour ly  wind  velocity  readings  taken  by  ship  PAPA  in  the  Gulf  of  Alaska 
over  a  15  year  period.  After  suitable  detrending  to  remove  1  year,  6  month, 
12  hour  and  6  hour  cyclic  trend  components,  a  first  order  autoregressive 
Gamma  model  (Lewis,  1981)  was  fitted  to  the  data,  the  use  of  this  model 
being  suggested  by  the  shape  of  the  (marginal)  histogram  of  the  data  (Figure 
6.1),  the  autocorrelation  function  (Table  6.1)  and  the  shape  of  the  log  of 
the  normalized  period ogram  of  the  data  (Figure  6.2).  After  detrending 
there  is  still  a  slight  6-hour  effect  ( p»21, 900 )  because  this  cycle  varies 
in  intensity  over  the  15  years;  in  what  follows  this  will  be  ignored  and  the 
data  will  be  treated  as  stationary. 

It  is  not  the  object  here  to  give  the  above  analysis  in  detail  but 
rather  to  give  an  alternative  analysis  of  the  data  using  NEAR  models;  this 
involves  a  preliminary  transformation  of  the  data  to  an  exponential 
marginal  distribution.  This  is  suggested  firstly  by  the  fact  that  a  Weibull 
distribution  is  commonly  used  by  meteorologists  for  wind  velocity  data  and 
secondly  by  the  fact  that  Weibull  and  Gamma  distributions  fit  the  ship  PAPA 
wind  velocity  data  equally  well  (Lewis  and  Hugus,  1982);  a  power 
transformation  of  the  Weibull  then  leads  immediately  to  the  desired 
exponential. 

This  transformation  is  preferred  to  the  more  usual  transformation  to 
normality  which,  as  we  shall  argue,  is  not  appropriate  in  this  case.  The 
data  is  in  fact  finely  discretized,  with  zero  values  being  inlcuded.  After 
a  power  transform  to  normality  the  zero  values  show  up  as  a  group  of  values 
3till  at  zero,  whereas  the  non-zero  values  are  shifted  away  from  zero  to 
form  the  normal  part  of  the  data  distribution.  This  zero  value  problem  is 
not  critical  with  the  power  transform  to  exponent. ia l ity  since  this 


distribution  gives  rise  to  a  high  proportion  of  zero  values. 

It  is  possible,  but  extremely  tedious,  to  smooth  out  the  discretization 
in  the  data.  Even  then,  however,  one  cam  anticipate  that  after  a 
trams format ion  to  marginal  normality  the  time  series  will  be  nonlinear; 
there  is  no  guarantee  that  such  a  transformation  will  produce  linearity. 
Thus  we  have  preferred  to  transform  for  marginal  exponent iality  and  attempt 
to  incorporate  nonlinearity  into  the  modelling. 

The  histogram  of  the  transformed  data  Xn'«Xn2,185,  is  shown  in  Figure 
6.3,  where  the  power  trams format ion  to  exponent iality  ham  been  determined 
iteratively  so  that  the  coefficient  of  variation  of  the  tramsformed  data  is 
unity.  This  transformation  does  affect  the  correlation  structure  of  the 
data,  as  shown  in  Table  6.1;  the  table  gives  comparisons  of  data  and  model 
autocorrelations  both  before  and  after  the  trams format ion. 


Figure  6.3.  Histogram  of  transformed,  detrended  ulnd  velocity  data.  The 
transformation  Is  Xn'«Xn*,x*s,  uhere  the  Xn,  n-1 ,2, ... ,43,800  are  the 
detrended  original  data.  The  shape  Is  clearly  exponential,  hut  the 
estimated  skeuness  and  kurtosls  have  values  greater  than  the  theoretical 
values  of  2.0  and  6.0,  caused  by  a  small  percentage  of  very  large  values. 


Table  6.1 


Pit  of  Wind  Speed  Date  Autocorrelation* 


Detrended  Series  Trams formed  Detrended  Series 


(Xn>  (Xn*-1®5) 


Pl,,# 

Est .  p(  t) 

1*1 

1*1 

Pi"1 

Est . p( 4 ) 

AR<2) 

1*1 

.8194 

^”-.8194 

.0000 

.0000 

.7902 

Pi"-. 7902 

.7902 

•  OOOO 

.6714 

P2"«.689B 

.0184 

.0346 

.6243 

P2"-“.  6589 

.6589 

.0000 

.5502 

P3”’.5635 

.0133 

.0391 

.4933 

P3“-.5324 

.5454 

.0130 

.4508 

P4"«. 4698 

.0190 

.0541 

.3898 

P4“-. 4439 

.4519 

.0080 

.3694 

V-.3764 

.0070 

.0431 

.3080 

P5"«.3511 

.3744 

.0233 

Table  6.1 

The  second 

column 

shous  the  estimated  autocorrelations 

for 

detrended 

series/  these  are 

close  in 

value  to 

the  powers 

of  Pj"-0 

.8194 

t he  first  column,  Indicating  a  flood  fit  (Column  3)  to  a  model  with  AR(l) 
autocorrelation  structure.  After  transformation,  this  A R(l)  fit  is  no 
longer  valid  (columns  4  and  5)  since  the  absolute  values  of  the  differences 
(Column  4)  are  nou  much  larger.  A  much  better  and  adequate  fit  is  obtained 
ult'i  A R(2)  autocorrelation  structure,  as  indicated  by  columns  6,  7  and  8. 


Column  2  in  Table  6.1  gives  the  estimated  autocorrelations,  p|"  of 
the  detrended  data;  the  standard  error  of  each  of  these  autocorrelation 
estimates  separately  is  given  approximately  by  l/(N)1/2«l/( 43, 800 l1/2 
“0.005.  The  first  column  of  Table  l  gives  the  fitted  autocorrelations  for  a 
model  with  AR(i)-type  autocorrelations,  just  pj."*  -  (0.8194)*,  for  lags 
f“l,2,...,5.  The  differences  (Column  3)  are  all  very  small  in  practical 
importance,  although  some  are  perhaps  statistically  significant  in  view  of 
the  large  sample  size. 

Column  6  in  Table  6.1  gives  the  estimated  autocorrelations,  p*“ ,  for 
the  transformed  data;  the  transformation  consistently  lowers  the 
autocorrelations.  However,  columns  5  and  7,  which  give  the  lilted  AR(  1 )  and 
AR(2)  correlation  values,  respectively,  show  that  a  model  with  the  AR(2) 
correlation  structure  is  definitely  preferable.  The  fit  is  borne  out  by  a 
periodogram  plot  (not  given),  and  the  analysis  will  be  continued  on  this 
basis. 

7bu8,  a  NEAR(2)  model  is  a  candidate  for  representing  the  transformed 
data,  and  if  p(l)  and  p(2)  are  fixed  at  the  estimated  values  of  *  0.7902 
and  p2”  “  0.6589,  then  the  corresponding  ai-a^Oi  and  a2«ot202  from  (5.4)  are, 
respectively,  0.7175  and  0.0920.  There  are  still  two  degrees  of  freedom 
left  in  fitting  the  model,  represented  by  choice  of  parameters  and  a2 
greater  than  or  equal  to  0.7175  and  0.0920  respectively,  with  cti+a2<l. 

Figure  6.4  shows  the  cumulative  periodogram  of  the  usual  AR(2)  model 
linear  residuals,  Rn“Xn’ -ajXn-i ' -a2Xn-2 ' ,  of  the  transformed  data.  This  is 
almost  straight  (ignoring  the  slight  effect  at  period  6  hours).  At  this 
point  it  might  be  thought  that  the  usual  second  order  autoregressive  model 
is  adequate.  We  shall  however  now  develop  an  extended  residual  analysis  for 
higher  order  dependence  which  explores  further  fitting  of  the  NEAR( 2 )  model 
to  the  transformed  data. 
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Figure 


p-(o/N  equals  cycles  per  3  hours) 
Cumulative  perlodogram  of  linear  residuals 


where 


Rn=Xn‘ -0.7J7SXn- j ’ -O . 0920Xn -2 ' ,  /or  the  detrended  and  transformed  wind 
velocity  data.  Since  the  cumulative  perlodogram  is  almost  straight,  the 
Indication  Is  that  the  residuals  are  uncorrelated .  Note  that  there  Is  still 
the  suggestion  of  a  six  hour  effect,  Indicated  by  the  final  upuard  turn  of 
the  graph . 


7 .  A  RESIDUAL  ANALYSIS  FOR  THE  NEAR(  2  )  MODEL 
7 . l  General  Results 

It  has  already  been  remarked  that  the  autocorrelations  p(  ( )  are 
insufficient  to  describe  the  dependency  structure  of  NEAR( 2 )  models.  A 
natural  next  step  might  be  to  examine  higher  order  joint  moments  and  their 
associated  spectra.  The  functions  so  obtained,  such  as  the  bispectrum, 
(Tukey  (1959)),  are  often  found  to  be  difficult  to  calculate  and  hard  to 
interpret.  Rather  than  follow  this  course,  a  residual  analysis  for 
nonlinear  autoregressive  models  is  proposed.  The  thrust  of  this  analysis 
is  that  the  standard  process  of  fitting  and  validating  a  linear 
autoregressive  model  should  be  carried  out  beyond  the  customary  final  stage 
at  which  uncorrelated  residuals  are  obtained  (as  in  the  previous  section). 
The  usual  presumption  is  that  the  residuals  are  not  only  uncorrelated  but 
also  independent.  This  need  not  be  the  case,  as  will  be  exemplified  for  the 
wind  velocity  series.  Also,  uncorrelated  but  dependent  residuals  are 
obtained  for  KEAR( 2 )  processes.  Thus  residuals  should  be  analysed  for 
further  dependency.  Any  found  is  then  evidence  that  a  standard  linear, 
constant  coefficient  second  order  autoregressive  model  is  deficient.  With 
norma] ly  distributed  time  series  data  this  might  suggest  that  Gaussian 
nonlinear  modelling  should  be  explored.  With  data  marginally  distributed 
in  some  other  identifiable  manner,  the  exploration  of  a  selected  type  of 
nonlinear  model  with  specified  marginal  distribution  and  autocorrelation 
function  is  suggested;  it  should  then  have  some  higher  order  dependency 
properties  of  its  autoregressive  residuals  in  agreement  with  those  of  the 
data.  This  is  the  course  exemplified  here.  The  proposed  residual  analysis 
is  further  explored  in  Lawrance  and  Lewis  (1984a). 

In  Section  9  higher  order  dependency  properties  of  the  uncorrelated 
residuals  from  the  wind  velocity  data  are  compared  with  those  derived  in 
Section  8  for  the  NEAR(2)  model.  This  stage  can  be  informative  from  both 


exploratory  and  estimation  considerations,  and  can  be  thought  o£  as  part  o£ 


the  model-refinement  process  common  to  much  statistical  methodology. 

It  might  be  thought  that  the  class  of  NEAR(  2 )  models  could  be  subjected 
to  a  residual  analysis  in  the  standard  manner  or  using  more  general  forms  of 
residuals  that  have  been  studied  by  Cox  and  Snell  (1968).  However, 
considering  the  near( 2 )  model  in  its  random  coefficient  form  (S.l),  the 
independent  innovation  is  now  trivariate,  consisting  of  (Kn,  Kn '  *  t*nEn  )  • 
No  way  of  estimating  this  trivaxiate  distribution  based  on  realized  (Xn) 
has  been  found,  even  assuming  knowledge  or  estimation  of  the  model 
parameters.  The  linear  autoregressive  residuals  are  never-the-less 
available,  and  given  as  in  the  previous  section  by 

RnmXn_alxn-l_a2xn-2 «  al~<*l0l«  a2“Q[2<32»  n  »  0,*1,±2, . . .  (7.1) 

We  now  show  that  these  residuals  are  uncorrelated  for  the  NEAR( 2 )  process. 
7.2  The  Residual  Theorem 

Theorem:  The  residuals  ( Rn,Rn+j )'  given  by  (7.1),  are  uncorrelated  for 

8»±1, *2 . 

Proof:  The  autocovariances  of  the  residuals  (7.1)  may,  for  8>1,  be 

written, 

Cov(Rn,  Rn+j )  =  Cov(Xn,  Rn+ 8 >~ alCov( Xn_ x ,  Rn+j )-a2Cov(  Xn_2 ,  Rn+8>  (7.2) 

-  Cov( Xn,  Rn+{ )-a^Cov( Xn,  Rn+4+1)-a2Cov(Xn,  Rn+e+2>' 
since  the  {Xn}  process  and  consequently  the  (Rn)  process  is  stationary.  The 
covariances  on  the  right  hand  side  are  all  of  the  same  type  and  given  by 
Cov(Xn,  Rti+  s  )  *  Cov(Xn,  <Xn+4  -  a^Xn+ -  a2Xn+4-2  ) } 

*  (Var(X)}(p(8)  -  aiP(8-l)  -  a2p(8-2)),  8-1,2,...  (7.3) 

An  identical  result  also  holds  for  8's  less  than  zero.  Hence  by  the  Yule- 
walker  equations  (5.4),  the  expression  in  brackets  is  zero,  and  so 

Corr(Rn,  Rn+I)  -  0,  8«±1,±2,...,  (7.4) 

as  was  to  be  proved.  That  these  residuals  are  uncorrelated  is  an  imnediate 


consequence  of  the  autocorrelations  following  Yule- Walker  equations;  this 
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emphasizes  that  residuals  of  the  fora  (7.1)  will  be  uncorrelated  for  any 
model  Whose  autocorrelations  satisfy  Yule-Walker  equations.  Equivalently, 
as  was  the  case  for  the  transformed  detrended  wind  velocity  data  (Figure 
6.4),  the  cumulative  periodogram  of  the  residuals  will  plot  linearly. 

Dependency  analysis  of  the  uncorrelated  residuals  (Rn),  n-1,2,...  could 
begin  with  scatter  plots  of  the  low-lag  adjacent  values;  any  patterns  or 
concentrations  are  suggestive  of  dependency.  Several  further  methods  for 
the  detection  of  dependency  could  be  proposed,  but  the  only  one  pursued  here 
involves  the  squares  of  residuals;  it  is  then  applied  in  Section  9  to  a 
continued  analysis  of  the  wind  velocity  data. 


8.  CROSSCOVARIANCE  ANALYSIS  OP  {Rn)  AND  (Bn*)  FOR  1BE  NEAR(2) 


After  the  satisfactory  fit  to  data  of  an  ordinary  linear  model,  the 
residuals,  Rn,  should  be  independent;  this  is  conveniently  investigated  by 
seeking  straight  cumulative  periodograms  for  the  residuals  and  for  the 
squared  residuals;  for  the  wind  velocity  data,  Figure  8.1  shows  that  the 
squared  residuals  have  an  obviously  curved  cumulative  periodogram.  Thus  a 
linear  AR(2)  model  for  this  data  is  definitely  not  adequate.  As  a  method 
for  probing  model  validity,  the  examination  of  squared  residuals  has 
previously  been  employed  by  McLeod  and  Li  (1983),  following  Granger  and 
Andersen  (1978);  these  latter  authors  investigated  bilinear  modelling  of 
uncorrelated  but  dependent  residuals  from  ARMA  models,  with  a  view  to 
improved  forecasting. 


Figure  8.1.  Cumulated  periodogram  values  for  the  squared  residuals,  Rnz, 
of  the  transformed  detrended  ulnd  velocity  data.  The  underlying  spectrum 
Is  clearly  not  flat/  the  Kolmogorov -Smirnov  test  statistic  Is  6.14  uhlch  Is 
very  high  compared  to  the  upper  .99  quantile  of  the  statistic  uhlch  Is 
1 .6 18. 
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Whilst  autocorrelations  of  the  squared  residuals  can  be  considered,  for 
the  NEAR( 2 )  model  this  involves  computation  of  36  terms,  mostly  distinct 
types  of  4th  order  moments.  A  more  tractable  suggestion  which  involves  only 
third  order  moments,  and  which  is  thus  the  next  step  up  after  ordinary 
autocorrelations,  is  to  use  the  croaacovartancaa  of  the  (Rn)  and  (Rn2) 
sequences;  apart  from  lag  o,  zero  values  would  be  found  for  linear  models. 
We  use  crosscovariances  rather  than  crosscorrelations  of  (Rn)  and  (Rn2) 
because  crosscorrelations  need  Var(Rn2)  which  involves  4th  order  moments  of 
the  xn  process. 

The  starting  point  for  the  calculation  of  Cov(Rn2,  Rn-|)  is  to  note  from 
the  definition  of  Rn  at  (7.1)  that 

E(Rn2Rn_4)  =  E(Rn2Xn_t)-a1E(Rn2Xn_|-1)-a2E(Rn2Xn_i-2),  1-0,  *1,  ±2,  .  .  .  (0.1) 

whence  there  is  the  structural  form, 

c°v( Rn2 ,  Rn-< )  *  Cov(Rn2,  Xn_j)-a1Cov(Rn2,  Xn_ t_ x )-a2Cov(  R„2 .  Xn_j_2 > •(  * 2 > 
Calculation  of  the  covariance  terms  in  (0.2)  requires  the  expanding  out  of 
Rn2,  taking  expectations  and  expression  in  terms  of  covariances.  Thus 
RnZxn- l  -  xnZxn-l  +  aizXn-l*xn-l  +  a2z*n-2z*n-l 

-  2aiXnXn_iXn_j  -  2a2XnXn_2xn- I 

+  2a^a2Xn— iXn_2Xn_ | ,  l-O, *1, *2, . . .  (0.3) 

The  conversion  to  covariance  yields 
Cov(  Rn2 ,  Xn_  | )  -  JX(«)  +  ai2Jx(  4-1)  +  a22J1(l-2) 

-  2ajj2(l)  -  2a2J3( t )  +  2aia2J2( t-1),  i«0,*l, *2, . . .  (0.4) 

Where 

JX(I)  -  Cov(Xn2,  Xn_ | ) ;  J2(  I )  -  Cov(XnXn_1,  Xn-f)> 

J3(l)  -  Cov(XnXn_2,  Xn_j).  (0.5) 

We  thus  see  the  types  of  third  order  joint  moments  for  the  (Xn)  process 
which  are  involved  in  the  Cov(Rn2,  Rn-t)  calculation;  these  are  general 
results . 

For  the  NEAR(  2 )  model  each  of  these  joint  moment s  has  to  be  obtained 


using  a  difference  equation,  considering  for  illustration  Ji(i)  for 


positive  lags,  square  each  side  of  the  NEAR(2)  defining  equation  (2.2)  and 
Multiply  by  Xn_f.  after  converting  to  the  required  covariances,  the 
recursion  is  found  to  be 

Jx(  *)  -  *101J1<  ,~1>  +  a2<»2Jl<  +  2(l-a1-a2)P(  «),  1-1,2,...  (8.6) 

This  equation  is  given  for  il lustration »  there  are  similar  equations  for 
J2(l)  and  J2(  I ),  and  various  special  cases.  The  complete  algorithm  for 
computing  Cov( Rn*,Rn_* )  for  I>1  for  the  NEAR( 2)  process  is  given  in  the 
Appendix. 

The  other  half  of  the  crosscovariance  function  of  (Rn}  and  {Rn2}  is 
Cov(Rn*,Rn_|)  for  I<1,  or  equivalently  Cov(Rn2,  Rn+i)  for  I>1>  we  now  show 
that  this  is  zero.  The  Key  result  in  establishing  this  fact  is  obtained  by 
first  defining  Gn_|  as  either  Xn-|2,  xn-lxn-l-l<  or  xn-lxn-l-2*  ^or  and 
showing  that  Cov(Rn,  G„_|)  for  I>1  are  all  zero. 

First  note  that 

cov(Rn.  Gn_j )  -  E(((xn-a1xn_i-a2xn-2)-*:(«n)HGn-|-E<Gn-l>)l 

-  EKXn-aiXn-i-azXn^HGn-t^Gn-i))).  (8-7) 

Now  substitute  for  Xn  from  (5.1)  to  obtain 

Cov(Rn,  Gn_()  -  E(((0iKn'-ai)Xn_i  +  ( 02xn"-a2  >^-2  +  *i»En)Ccn-l"'Etcn-l  ^  1  • 

(8.8) 

Since  ( Kn’  ,Kn" )  are  independent  of  ( Xn_  i,Xn_2 )  and  I^En  is  independent  of 
Gn_|,  the  right  hand  side  of  (8.8)  may  be  written 

Cov(Rn,  Gn-f)  -  E(/31Kn’-a1  )Cov(Xn_i,  Gn_j)+  E(  02xn"~a2  xn-2 »  Gn-I >• 

(8.9) 

By  the  definition  of  (Kn',Kn")  at  (5.3),  E(0xxn  -al)  an<*  E(02xn  -a2)  are 
both  zero,  and  hence 

Cov(  Rfj,  Gn_  | )  -  0 ,  1-1,2,...  (8.10) 

Finally,  we  note  that 
Cov(Rn2,  Rn+| )  -  Cov(Rn,  Rn-!2) 

-  Cov(R„,  Xn_*2)  +  a!2Cov(Rn,  Xn-i-x2)  +  a2*Cov(Rn,  Xn_f_x2) 

-  2ajCov(Rn,  Xn_ |Xn_ j_ x )  —  2a2cov(Rn,  xn- lxn- 1-2 ) 

+  2axa2Cov(Rn,  Xn- 1- ixn- 1-  2 ) • 


(8.11) 


By  (S.io)  all  the  covariances  in  (8.11)  are  sero,  and  the  desired  result  is 
estahl ished . 

in  a  similar  way  it  may  be  seen  that  Cov(Rnr,Rlv(.| )  for  I>1  and  all 
positive  integers  r  are  also  zeroj  this  does  not,  however,  imply  that  Ro  and 
Rn+I  are  independent.  All  joint  moments  would  have  to  be  sero  but  in 
particular  this  is  not  the  case  for  Cov(Rnz,  Rn-|)>  1*1,  as  results  in  the 
Appendix  indicate. 

It  should  perhaps  also  be  noted  explicitly  that  in  using  the  residuals 
{R„}  that  the  coefficients  a^  and  as  will  have  been  estimated;  this  may  in 
fact  induce  some  correlation  between  Rn2  and  Rn+f,  but  with  long  series  of 
data  the  effect  should  be  very  small. 


9. 


ANALYSIS  AND  MODELLING  OP 


WIND  VELOCITY  DATA 


Dependence  in  the  uncorrelated  second  order  residuals  of  the 
transformed  detrended  wind  velocity  data  has  already  been  demonstrated  by 
Figure  8.1.  Further  evidence  of  this  is  provided  by  the  non-sero 
crosscovariances  of  ( R,,,  Rn+i2)  given  in  Figure  9.1.  The  corresponding 
theoretical  crosscovariances  for  the  NEAR(2)  model  will  next  be  presented, 
having  been  computed  using  the  algorithm  described  in  Section  8  and  the 
Appendix. 


Crosscovorionces  leg 


Figure  9.1.  Croaacovartancea  between  the  second  order  reaiduala,  Rn,  and 
the  squared  aecond  order  reaiduala,  Rn+ j1,  for  the  trana formed  detrended 
utnd  velocity  data.  Sample  alee  la  N~43 , BOO .  There  la  a  atrong  negative 
value  at  lag  minus  one  and  pronounced  poaltlve  values  at  the  flrat  feu 
poattlve  1 age. 


At  this  point,  it  will  be  recalled,  the  NFAR(  2 )  model  has  not  been 
fitted  in  terms  of  all  4  parameters)  the  residuals  involve  the  model 
parameters  only  through  aj.  -  <>idl  and  a2  "  a202  and  <x\,0ii<x2'^2  have  not 
been  separately  estimated.  In  the  present  rather  exploratory  analysis  the 
estimation  problem  will  be  circumvented;  the  crosscovariances  of  (Rn,Rn+|*) 
will  be  given  for  four  representative  sets  of  parameter  values  in  the 
reduced  allowable  region,  as  constrained  by  (<*i»*i*a2>a2*al+a2<1 ), 
according  to  (5.6).  For  the  transformed  detrended  wind  velocity  data, 
ai~0.717S  and  a2=0.0920,  and  the  four  chosen  sets  of  ai,a2,0i,02  together 
with  their  associated  (p2,p3)  and  (b2,b3)  sets  are  given  in  Table  9.1. 

Table  9.1 

Parameters  Sets  for  NEAR(2)  Calculations 


<*1 

a2 

01 

02 

P2 

P3 

b2 

*>3 

A 

.72 

.0920 

.9965 

1.0000 

.0038 

.9962 

.9996 

.1874 

B 

.76 

.0920 

.9441 

1.0000 

.0585 

.9415 

.9940 

.1406 

C 

.77 

.1400 

.9318 

0.6571 

.0542 

.8610 

.7008 

.0786 

D 

.88 

.0920 

.8153 

1.0000 

.1744 

.8256 

.9826 

.0232 

In  Figure  9.2  the  crosscovariances  of  (Rn,  Rn+t*)  for  each  of  these  four 
cases  are  presented;  the  NEAR(2)  process  has  been  taken  with  mean  of  1.31  so 
as  to  correspond  to  the  transformed  detrended  wind  velocity  data.  First  it 
will  be  recalled  that  for  the  usual  linear  AR(2)  process,  these 
crosscovariances  would  all  be  zero,  apart  from  that  one  at  lag  zero  Which  is 
giving  a  non-standard  measure  of  skewness  of  the  residuals;  Case  D  is 
nearest  to  this  situation,  although  the  model  cannot  reduce  to  the  linear 
AR(  2 )  model;  it  can  reduce  to  the  linear  AR(1)  model  and  case  D  is  also 
nearest  to  this  situation. 

Case  C  is  nearest  in  qualitative  behaviour  to  the  crosscovariances  of 
the  wind  velocity  data  given  in  Figure  9.1,  but  as  we  have  seen,  cannot 
exhibit  the  negative  crosscovariance  shown  at  lag  minus  one.  This  is  an 


important  deficiency  of  modelling  by  the  NEAR(  2 )  process  although  matters 
would  have  been  considerably  worse  for  data  with  a  double  sided  CovCR,,2, 
Rn+l)  function.  The  other  cases  (A  and  B)  illustrate  some  of  the  diversity 
of  behaviour  producible  by  the  NEAR(  2 )  model)  more  investigation  of  these 
apaects  would  be  valuable.  Further  analysis  would  require  formal 
estimation  of  all  four  parameters,  or  perhaps  formal  estimation  of  (a^.a^) 
after  the  fixing  of  aj  and  &2  at  their  values  determined  by  estimates  of 
p(  1 )  and  p(  2  ) . 
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Figure  9.2.  Crosscovariances ,  for  the  NRA R(2)  model,  between  the  second 
order  residuals ,  Rn,  and  second  order  residuals  Rn+f2,  tilth  aj^O .7175  and 
02*0.0920,  The  sets  of  parameter  values  (,ai»<*2'01'02'P2'P3'*>2'*>3^  ar 0 
described  in  Table  9.1;  case  C  is  closest  to  that  of  the  uind  velocity  data 
given  tn  Figure  9.1.  Case  D  is  nearest  to  the  linear  EAR(l)  model.  The 
other  too  cases  Illustrate  more  of  the  diversity  of  behaviour  which  can  be 
shoun  by  this  NEAR(2)  crosscovariance  function.  The  lag  sero  covariance  is 
a  measure  (not  the  usual  one)  of  skeuness  of  the  residuals. 


10 .  OQMCLBSIGMS  AMD  FURTHER  ANALYSIS 

The  very  broad  four  parameter  NEAR( 2 )  time  series  model  having 
exponential  marginals  and  the  correlation  structure  of  a  linear  AR( 2 )  model 
has  been  established.  A  preliminary  fit  of  the  NEAR(2)  model  has  been  made 
to  a  very  long  series  of  wind  speed  data,  the  data  having  been  detrended  and 
transformed  so  as  to  have  exponentially  distributed  marginals.  A  residual 
analysis  has  been  based  on  the  crosscovariances  between  the  residuals  and 
squared  residuals,  and  its  utility  in  probing  higher  order  dependence  in 
the  ( Xn }  process  has  been  demonstrated;  in  particular  it  has  highlighted  a 
strong  difference  in  higher  order  crosscovariance  of  the  data  and  the 
NEAR( 2 )  model  at  lag  minus  one. 

A  possible  extension  of  the  NRAR(2)  model  which  retains  the  marginal 
exponential  distribution  of  the  data  is  obtained  by  noting  that  the  theorem 
in  Section  3  does  not  require  independence  of  the  random  coefficient 
sequences  {Kn' }  and  {Kn" }.  By  allowing  these  to  be,  say,  moving  average 
sequences,  the  higher  order  structure  of  the  {Xn}  process  might  be  extended 
to  accomodate  the  negative -valued  spike  in  Figure  9.1  at  lag  minus  one. 
This  approach  has  been  taken  by  McKenzie  ( 1981 )  but  it  is  not  known  how 
tractable  the  resulting  structure  of  the  NEAR( 2 )  model  would  be. 

Finally,  we  remark  that  a  likelihood  conditional  on  the  first  two 
values,  Xi  and  X2,  can  be  written  down  for  the  process.  However  the 
likelihood  function  is  difficult  to  use  because  it  becomes  infinite  at 
parameter  values  on  the  border  of  the  parameter  space  corresponding  to  an 
EAR(  ? )  model.  These  considerations  suggest  that  a  reduced  three  parameter 
model  is  needed  which  avoids  these  singularities  but  still  allows  for  a 
higher  order  dependency  in  the  model.  Such  a  likelihood  analysis  would  also 
need  model  validation  which  could  be  based  on  the  higher  order  residual 
analysis  presented  here. 

An  extension  of  the  residual  analysis  given  here  based  on  reversed 


residuals  is  possible  (Lawrance  and  Lewis,  1984b);  however,  these  were  not 
used  because  residual  analysis  has  already  turned  up  discrepancies  between 


the  data  and  the  model  which  need  further  explanation. 


appendix 

Algorithm  for  Commuting  the  non  negative  half  off  the  (Residual.  Residual 
Squared!  Croaeoovarianoe  Function  for  the  WEAR/ 2 )  Model 

Input  <*i ,  a2 , 03 , 02  *  L 

0.  ai~<xi0i/  a2-<X202  > 

p(0)  *  l;  p(l)  *  ai/(l-a2);  p(2)  -  aip(l)  4  a2 

p(<)  «  aip(l-l)  4  a2p(l-2),  for  1-2, 3,..., L. 
l.  Jx(0)  *  4;  J!(1)  *  [4a1(014a202 )  +  2( l-a!-a2 ) )/( i-a2*02 ) . 

2a.  Ji(«)  -  aiJ1(*41)  4  a2Ji(  4+2)  for  1—1, -2. 

2b.  J\(i)  -  a101J1(l-l)  4  a202Jl(i-2)  +  2( l-ai-a2 )P( I )  for  <-2,3,...,L. 

3.  J2( 0 )  -  Ji(l)  -  p(l)  +  1;  J2(l)  *  Ji(-l)  -  p<l)  +  1. 

4a.  J2(-i )  *  aiJ2(0)  +  a2J2(l). 

4b.  J2(<)  =  aiJx(«-l)  +  a2J2<*-l)  +  2aL  +  [l+p<l)]a2 

-[l+p(l)]  +  [1  +  p(  l-l)](l-ai-a2)  for  <-2,3,...,L. 

5.  J3(0)  -  Ji(2)  -  p(2)  +  1;  J3<1)  -  J2(2)+p(l)  -  p(2); 

J3<2)  =  Ji(  -2)  -  p(  2)  +  1. 

6.  J3(<)  -  aiJ2(«-i)  +  a2Ji(«-2)  +  [1+P(i)]ai  +  2a2 

+  [l+p(  f-2)](l-a1-a2)  -  [ l+p(  2)]  for  <-3,4,5, ...  ,L. 

7.  J(<)  -  Jid)  +  ax*Ji(  1-1)  +  a22Ji(<-2)  -  2aiJ2(<) 

-  2a2J3( t )  +  2aia2J2( 1-1 )  for  <-o,l,...,L. 

L. 


8.  Cov( ,  Rn-1)  -  J(l)  -  aiJ(l+l)  -  a2J( <+2 )  for  1-0,1,..., 
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